### This script creates the analysis and figures in "Substantive Divergence: 
### The Meaning of Public Opinion on Government Spending in Red and Blue"
### Katherine Krimmel and Kelly Rader
### Aug 14, 2020

### NOTE: Figures are saved to the working directory as pdfs.
### Tables are saved to the working directory as tex files.

library(car)
library(MASS)
library(sandwich)
library(lmtest)
library(ggplot2)
library(ggthemes)
library(vcd)
library(plyr)
library(rms)
library(nnet)
library(apsrtable)
library(aod)
library(multcomp)
library(pastecs)
library(foreign)
library(reshape2)
library(grDevices)
library(gridExtra)


survey <- read.csv("substantive_divergence_replication_data.csv", header=TRUE, sep=",",stringsAsFactors=F)


##################################### 
##################################### mean responses for plots
##################################### 
overallmeans <-  cbind.data.frame(sort(100*c(apply(survey[grep("(?=.*sO_)",names(survey),value=T,perl=T)],2,mean,na.rm=T))),
                                  c("Taxes","Research","Social Security and Medicare","Businesses","Immigrants", "Other domestic",
                                    "Debt","Infrastructure","General domestic",
                                    "Education","Healthcare","Government employees","Waste","Social services",
                                    "Means-tested programs","Foreign aid","Defense","General"))

overallmeans$sO <- rownames(overallmeans)
names(overallmeans)<-c("allpct","name","cat")
overallmeans$mega <- NA
overallmeans$mega[overallmeans$cat=="sO_socserv"|overallmeans$cat=="sO_welf"|overallmeans$cat=="sO_health"|overallmeans$cat=="sO_educ"|
                    overallmeans$cat=="sO_otherdom"|overallmeans$cat=="sO_domestic"|overallmeans$cat=="sO_infr"|overallmeans$cat=="sO_immi"|
                    overallmeans$cat=="sO_biz"|overallmeans$cat=="sO_ss"|overallmeans$cat=="sO_rd"] <- "Domestic"
overallmeans$mega[overallmeans$cat=="sO_defense"|overallmeans$cat=="sO_foraid"] <- "Foreign"
overallmeans$mega[overallmeans$cat=="sO_every"] <- ""
overallmeans$mega[overallmeans$cat=="sO_waste"|overallmeans$cat=="sO_govemp"|overallmeans$cat=="sO_debt"|overallmeans$cat=="sO_tax"] <- "Other"


##################################### 
##################################### FIGURE 1: MAIN OPEN-ENDED ANSWER BARPLOT
##################################### 

plotdata <- overallmeans
names(plotdata)[names(plotdata)=="allpct"] <- "pct"

pdf("main_open_ended.pdf",height=8,width=9.5)
ggplot(plotdata, aes(x=reorder(name,pct), fill=mega, y=pct)) + geom_bar(stat="identity") +
  scale_fill_grey(guide=FALSE, start=.4) +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  geom_text(aes(label=paste(round(pct,digits=0),"%",sep="")), hjust=1.2, size=5) + theme_bw() + coord_flip() +
  scale_y_discrete(breaks=NULL) +
  theme(plot.title=element_text(size=13, vjust=2, hjust=.5), panel.grid.major = element_blank(),
        axis.title.x = element_blank(), axis.title.y = element_blank(), panel.border=element_blank(), 
        axis.text.y = element_text(size=15), axis.ticks = element_blank(), strip.background = element_blank(),
        strip.text = element_text(size=15, face="bold"))
dev.off()  


##################################### 
##################################### FIGURE 2: SCATTERPLOT COMPARING DEMOCRATS AND REPUBLICANS
##################################### 

survey_rep <- subset(survey,rep==1)
plotdata_r <- cbind.data.frame(100*c(apply(survey_rep[grep("(?=.*sO_)",names(survey_rep),value=T,perl=T)],2,mean,na.rm=T)))
plotdata_r$sO <- rownames(plotdata_r)                            
names(plotdata_r)<-c("pct_r","cat")
plotdata_r <- merge(plotdata_r,overallmeans, by="cat")
plotdata_r <- plotdata_r[order(plotdata_r$allpct),]

survey_dem <- subset(survey,dem==1)
plotdata_d <- cbind.data.frame(100*c(apply(survey_dem[grep("(?=.*sO_)",names(survey_dem),value=T,perl=T)],2,mean,na.rm=T)))
plotdata_d$sO <- rownames(plotdata_d)                            
names(plotdata_d)<-c("pct_d","cat")

plotdata_scatter <- merge(plotdata_d, plotdata_r, by="cat")
plotdata_scatter <- subset(plotdata_scatter, name!="General")

pdf("scatter_open_ended.pdf",height=8,width=8)
ggplot(plotdata_scatter, aes(x=pct_d, y=pct_r)) +
 geom_point()+
  geom_text(aes(x=pct_d+0.25, y=pct_r-0.25, label=name), size=3, hjust=0)+
  geom_abline(intercept=0, slope=1,linetype="dashed")+
  xlim(0, 31) + ylim(0, 31)+
  xlab("\n Democrats \n") + ylab("\n Republicans \n")+
  ggtitle("\n Percentage of Democrats and Republicans Mentioning Each Type of Spending \n")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5))+
  guides(size=FALSE)
dev.off()  

##################################### 
##################################### ASSOCIATION RESPONSE MODELS w/CLEAVAGES
##################################### 

dv <- c("sO_every", "sO_defense", "sO_foraid", "sO_welf", "sO_waste", "sO_socserv", "sO_health", "sO_educ", "sO_govemp", "sO_domestic", "sO_debt", "sO_infr", "sO_immi", "sO_otherdom", "sO_biz", "sO_ss", "sO_rd", "sO_tax")

# to store party results, lpm
party_coefs <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs <- merge(party_coefs,overallmeans,by.x="dv",by.y="cat")
# to store party results, logit
party_coefs_l <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_l) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_l <- merge(party_coefs_l,overallmeans,by.x="dv",by.y="cat")

for(i in 1:length(dv)){
  # party simple model, lpm w/HC3 ses
  model <- paste("model",party_coefs$dv[i],"_rep",sep="")
  m <- lm(as.formula(paste(party_coefs$dv[i],"~ rep + ind")), data=survey)
  assign(model,m)
  party_coefs[i,2] <- coef(eval(parse(text=(paste("model",party_coefs$dv[i],"_rep", sep="")))))[2]
  party_coefs[i,3] <- sqrt(diag(vcovHC(eval(parse(text=(paste("model",party_coefs$dv[i],"_rep", sep="")))))))[2]
  # party simple model, logit
  model <- paste("model",party_coefs_l$dv[i],"_rep_l",sep="")
  m <- glm(as.formula(paste(party_coefs_l$dv[i],"~ rep + ind")), data=survey, family=binomial)
  assign(model,m)
  party_coefs_l[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_l$dv[i],"_rep_l", sep="")))))[2]
  party_coefs_l[i,3] <- summary(eval(parse(text=(paste("model",party_coefs_l$dv[i],"_rep_l", sep="")))))$coefficients[2,2]
  # all model, lpm w/HC3 ses
  modelmv <- paste("modelmv",party_coefs$dv[i], sep="")
  mv <- lm(as.formula(paste(party_coefs$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey)
  assign(modelmv,mv)
  party_coefs[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs$dv[i], sep="")))))[2]
  party_coefs[i,5] <- sqrt(diag(vcovHC(eval(parse(text=(paste("modelmv",party_coefs$dv[i], sep="")))))))[2]
  # all model logit
  modelmv_l <- paste("modelmv",party_coefs_l$dv[i],"_l",sep="")
  mv <- glm(as.formula(paste(party_coefs_l$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, family=binomial)
  assign(modelmv_l,mv)
  party_coefs_l[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_l$dv[i],"_l", sep="")))))[2]
  party_coefs_l[i,5] <- summary(eval(parse(text=(paste("modelmv",party_coefs_l$dv[i],"_l", sep="")))))$coefficients[2,2]
  }

##################################### 
##################################### FIGURE 3: COEF PLOT
##################################### 

pdf("party_comparison.pdf",height=7,width=6)
ggplot(party_coefs, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                           ymax = repcoef + qnorm(1 - 0.05 / 2) * repse))+
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
     
  geom_pointrange(data=party_coefs,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                          ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents") +
  ylim(-.22,.22) +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()



##################################### 
##################################### ASSOCIATION RESPONSE MODELS, PARTY by SP_MAIN ANSWER
##################################### 

dv <- c("sO_every", "sO_defense", "sO_foraid", "sO_welf", "sO_waste", "sO_socserv", "sO_health", "sO_educ", "sO_govemp", "sO_domestic", "sO_debt", "sO_infr", "sO_immi", "sO_otherdom", "sO_biz", "sO_ss", "sO_rd", "sO_tax")

###### to store too much, lpm results
party_coefs_toomuch <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_toomuch) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_toomuch <- merge(party_coefs_toomuch,overallmeans,by.x="dv",by.y="cat")

###### to store too much, logit results
party_coefs_toomuch_l <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_toomuch_l) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_toomuch_l <- merge(party_coefs_toomuch_l,overallmeans,by.x="dv",by.y="cat")

for(i in 1:length(dv)){
  # party simple model, lpm
  model <- paste("model",party_coefs_toomuch$dv[i],"_rep",sep="")
  m <- lm(as.formula(paste(party_coefs_toomuch$dv[i],"~ rep + ind")), data=survey, subset=sp_main==3)
  assign(model,m)
  party_coefs_toomuch[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_toomuch$dv[i],"_rep", sep="")))))[2]
  party_coefs_toomuch[i,3] <- sqrt(diag(vcovHC(eval(parse(text=(paste("model",party_coefs_toomuch$dv[i],"_rep", sep="")))))))[2]
  # party simple model, logit
  model <- paste("model",party_coefs_toomuch_l$dv[i],"_rep_l",sep="")
  m <- glm(as.formula(paste(party_coefs_toomuch_l$dv[i],"~ rep + ind")), data=survey, subset=sp_main==3, family=binomial)
  assign(model,m)
  party_coefs_toomuch_l[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_toomuch_l$dv[i],"_rep_l", sep="")))))[2]
  party_coefs_toomuch_l[i,3] <- summary(eval(parse(text=(paste("model",party_coefs_toomuch_l$dv[i],"_rep_l", sep="")))))$coefficients[2,2]
  # all model, lpm
  modelmv <- paste("modelmv",party_coefs_toomuch$dv[i], sep="")
  mv <- lm(as.formula(paste(party_coefs_toomuch$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==3)
  assign(modelmv,mv)
  party_coefs_toomuch[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_toomuch$dv[i], sep="")))))[2]
  party_coefs_toomuch[i,5] <- sqrt(diag(vcovHC(eval(parse(text=(paste("modelmv",party_coefs_toomuch$dv[i], sep="")))))))[2]
  # all model, logit
  modelmv <- paste("modelmv",party_coefs_toomuch_l$dv[i],"_l",sep="")
  mv <- glm(as.formula(paste(party_coefs_toomuch_l$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==3, family=binomial)
  assign(modelmv,mv)
  party_coefs_toomuch_l[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_toomuch_l$dv[i],"_l", sep="")))))[2]
  party_coefs_toomuch_l[i,5] <- summary(eval(parse(text=(paste("modelmv",party_coefs_toomuch_l$dv[i],"_l", sep="")))))$coefficients[2,2]
}


###### too little
party_coefs_toolittle <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_toolittle) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_toolittle <- merge(party_coefs_toolittle,overallmeans,by.x="dv",by.y="cat")

###### too little logit
party_coefs_toolittle_l <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_toolittle_l) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_toolittle_l <- merge(party_coefs_toolittle_l,overallmeans,by.x="dv",by.y="cat")


for(i in 1:length(dv)){
  # party simple model
  model <- paste("model",party_coefs_toolittle$dv[i],"_rep",sep="")
  m <- lm(as.formula(paste(party_coefs_toolittle$dv[i],"~ rep + ind")), data=survey, subset=sp_main==1)
  assign(model,m)
  party_coefs_toolittle[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_toolittle$dv[i],"_rep", sep="")))))[2]
  party_coefs_toolittle[i,3] <- sqrt(diag(vcovHC(eval(parse(text=(paste("model",party_coefs_toolittle$dv[i],"_rep", sep="")))))))[2]
  # party simple model, logit
  model <- paste("model",party_coefs_toolittle_l$dv[i],"_rep_l",sep="")
  m <- glm(as.formula(paste(party_coefs_toolittle_l$dv[i],"~ rep + ind")), data=survey, subset=sp_main==1, family=binomial)
  assign(model,m)
  party_coefs_toolittle_l[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_toolittle$dv[i],"_rep_l", sep="")))))[2]
  party_coefs_toolittle_l[i,3] <- summary(eval(parse(text=(paste("model",party_coefs_toolittle$dv[i],"_rep_l", sep="")))))$coefficients[2,2]
  # all model
  modelmv <- paste("modelmv",party_coefs_toolittle$dv[i], sep="")
  mv <- lm(as.formula(paste(party_coefs_toolittle$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==1)
  assign(modelmv,mv)
  party_coefs_toolittle[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_toolittle$dv[i], sep="")))))[2]
  party_coefs_toolittle[i,5] <- sqrt(diag(vcovHC(eval(parse(text=(paste("modelmv",party_coefs_toolittle$dv[i], sep="")))))))[2]
  # all model, logit
  modelmv <- paste("modelmv",party_coefs_toolittle_l$dv[i],"_l",sep="")
  mv <- glm(as.formula(paste(party_coefs_toolittle_l$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==1, family=binomial)
  assign(modelmv,mv)
  party_coefs_toolittle_l[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_toolittle_l$dv[i],"_l", sep="")))))[2]
  party_coefs_toolittle_l[i,5] <- summary(eval(parse(text=(paste("modelmv",party_coefs_toolittle_l$dv[i],"_l", sep="")))))$coefficients[2,2]
  
}


##################################### 
##################################### FIGURE 4: COEF PLOT by SP_MAIN RESPONSE
#####################################

###### too much
party_toomuch_graph <- ggplot(party_coefs_toomuch, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                        ymax = repcoef + qnorm(1 - 0.05 / 2) * repse))+
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
  
  geom_pointrange(data=party_coefs_toomuch,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                       ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents\nToo Much") +
  ylim(-.33,.33) +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_blank(),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))


###### too much
party_toolittle_graph <- ggplot(party_coefs_toolittle, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                                    ymax = repcoef + qnorm(1 - 0.05 / 2) * repse)) +
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
  
  geom_pointrange(data=party_coefs_toolittle,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                                   ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents\nToo Little") +
  ylim(-.33,.33) +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_blank(), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))


############ Putting the too much and too little graphs together
pdf("party_toomuch_toolittle.pdf",height=9,width=10)
grid.newpage()
grid.draw(cbind(ggplotGrob(party_toomuch_graph), ggplotGrob(party_toolittle_graph), size = "last"))
dev.off()


##################################### 
##################################### 
##################################### APPENDIX
#####################################
##################################### 


######################################## 
######################################## APPENDIX: FIG A2, weighted vs unweighted
######################################## 

## mean mentioned, weighted
overallmeans_w <-  
  cbind.data.frame(sort(100*c(apply(survey[grep("(?=.*sO_)",names(survey),value=T,perl=T)],2,
                                    function(x) weighted.mean(x,survey$wgt,na.rm=T)))))

overallmeans_both <- merge(overallmeans_w, overallmeans, by="row.names")


names(overallmeans_both)<-c("cat","allpct_w","allpct","name")
overallmeans_both$mega <- NA
overallmeans_both$mega[overallmeans_both$cat=="sO_socserv"|overallmeans_both$cat=="sO_welf"|overallmeans_both$cat=="sO_health"|overallmeans_both$cat=="sO_educ"|
                         overallmeans_both$cat=="sO_otherdom"|overallmeans_both$cat=="sO_domestic"|overallmeans_both$cat=="sO_infr"|overallmeans_both$cat=="sO_immi"|
                         overallmeans_both$cat=="sO_biz"|overallmeans_both$cat=="sO_ss"|overallmeans_both$cat=="sO_rd"] <- "Domestic"
overallmeans_both$mega[overallmeans_both$cat=="sO_defense"|overallmeans_both$cat=="sO_foraid"] <- "Foreign"
overallmeans_both$mega[overallmeans_both$cat=="sO_every"] <- ""
overallmeans_both$mega[overallmeans_both$cat=="sO_waste"|overallmeans_both$cat=="sO_govemp"|overallmeans_both$cat=="sO_debt"|overallmeans_both$cat=="sO_tax"] <- "Other"

##################################### MAIN OPEN-ENDED ANSWER BARPLOT, weighted and unweighted

plotdata <- cbind.data.frame(c(overallmeans_both$allpct,overallmeans_both$allpct_w),
                             c(rep("un",18),rep("wgt",18)),
                             c(as.character(overallmeans_both$name),as.character(overallmeans_both$name)),
                             c(as.character(overallmeans_both$mega),as.character(overallmeans_both$mega)))

names(plotdata) <- c("pct","uw","name","mega")

pdf("main_open_ended_wgt_comparison.pdf",height=8,width=9.5)
ggplot(plotdata, aes(x=reorder(name,pct), y=pct, fill=uw)) + geom_bar(stat="identity", position="dodge") +
  scale_fill_grey(guide=FALSE, start=.4) +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  theme_bw() + coord_flip() +
  scale_y_continuous(sec.axis = dup_axis()) +
  theme(plot.title=element_text(size=13, vjust=2, hjust=.5), panel.grid.major = element_blank(),
        axis.title.x = element_blank(), axis.title.y = element_blank(), panel.border=element_blank(), 
        axis.text.y = element_text(size=15), axis.ticks.y = element_blank(), strip.background = element_blank(),
        strip.text = element_text(size=15, face="bold"))
dev.off() 

### note, weighted is light gray, unweighted is dark gray



######################################## 
######################################## APPENDIX: "right amount" models 
######################################## 

###### to store right amount, lpm results
party_coefs_right <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_right) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_right <- merge(party_coefs_right,overallmeans,by.x="dv",by.y="cat")

###### to store too much, logit results
party_coefs_right_l <- cbind.data.frame(dv,rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)),rep(NA,length(dv)))
names(party_coefs_right_l) <- c("dv","repcoef","repse","repcoefmv","repsemv")
party_coefs_right_l <- merge(party_coefs_right_l,overallmeans,by.x="dv",by.y="cat")

for(i in 1:length(dv)){
  # party simple model, lpm
  model <- paste("model",party_coefs_right$dv[i],"_rep",sep="")
  m <- lm(as.formula(paste(party_coefs_right$dv[i],"~ rep + ind")), data=survey, subset=sp_main==2)
  assign(model,m)
  party_coefs_right[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_right$dv[i],"_rep", sep="")))))[2]
  party_coefs_right[i,3] <- sqrt(diag(vcovHC(eval(parse(text=(paste("model",party_coefs_right$dv[i],"_rep", sep="")))))))[2]
  # party simple model, logit
  model <- paste("model",party_coefs_right_l$dv[i],"_rep_l",sep="")
  m <- glm(as.formula(paste(party_coefs_right_l$dv[i],"~ rep + ind")), data=survey, subset=sp_main==2, family=binomial)
  assign(model,m)
  party_coefs_right_l[i,2] <- coef(eval(parse(text=(paste("model",party_coefs_right_l$dv[i],"_rep_l", sep="")))))[2]
  party_coefs_right_l[i,3] <- summary(eval(parse(text=(paste("model",party_coefs_right_l$dv[i],"_rep_l", sep="")))))$coefficients[2,2]
  # all model, lpm
  modelmv <- paste("modelmv",party_coefs_right$dv[i], sep="")
  mv <- lm(as.formula(paste(party_coefs_right$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==2)
  assign(modelmv,mv)
  party_coefs_right[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_right$dv[i], sep="")))))[2]
  party_coefs_right[i,5] <- sqrt(diag(vcovHC(eval(parse(text=(paste("modelmv",party_coefs_right$dv[i], sep="")))))))[2]
  # all model, logit
  modelmv <- paste("modelmv",party_coefs_right_l$dv[i],"_l",sep="")
  mv <- glm(as.formula(paste(party_coefs_right_l$dv[i],"~ rep + ind + scale(age) + income5 + income2 + income3 + income4 + black + latino + otherrace + scale(edu) + male")), data=survey, subset=sp_main==2, family=binomial)
  assign(modelmv,mv)
  party_coefs_right_l[i,4] <- coef(eval(parse(text=(paste("modelmv",party_coefs_right_l$dv[i],"_l", sep="")))))[2]
  party_coefs_right_l[i,5] <- summary(eval(parse(text=(paste("modelmv",party_coefs_right_l$dv[i],"_l", sep="")))))$coefficients[2,2]
}



######################################## 
######################################## APPENDIX: FIG A3, logit coeff plot
######################################## 

pdf("party_comparison_logit.pdf",height=7,width=6)
ggplot(party_coefs_l, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                          ymax = repcoef + qnorm(1 - 0.05 / 2) * repse))+
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
  
  geom_pointrange(data=party_coefs_l,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                         ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents, Logit") +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()


######################################## 
######################################## APPENDIX: FIG A4, "too much" logit coeff plot
######################################## 


pdf("party_comparison_toomuch_logit.pdf",height=7,width=6)
ggplot(party_coefs_toomuch_l, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                                  ymax = repcoef + qnorm(1 - 0.05 / 2) * repse))+
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
  
  geom_pointrange(data=party_coefs_toomuch_l,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                                 ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents, Logit\nToo Much") +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()




######################################## 
######################################## APPENDIX: FIG A5, "right amount" coeff plot 
######################################## 

pdf("party_comparison_rightamount.pdf",height=7,width=6)
ggplot(party_coefs_right, aes(x=reorder(name, allpct), y=repcoef, ymin = repcoef - qnorm(1 - 0.05 / 2) * repse,
                              ymax = repcoef + qnorm(1 - 0.05 / 2) * repse))+
  geom_pointrange(size=.6, position=position_nudge(x = .1)) +
  
  geom_pointrange(data=party_coefs_right,aes(x=name, y=repcoefmv, ymin = repcoefmv - qnorm(1 - 0.05 / 2) * repsemv,
                                             ymax = repcoefmv + qnorm(1 - 0.05 / 2) * repsemv), 
                  size=.6, position=position_nudge(x = -.1), color="gray75") +
  facet_grid(mega ~ ., scales="free_y", space="free_y") +
  ggtitle("Republican vs. Democratic Respondents\nRight Amount") +
  ylim(-.33,.33) +
  geom_hline(yintercept = 0, lwd = .75, lty="dotted", alpha = .4) +
  coord_flip() + theme_bw() +
  theme(plot.title = element_text(size = 13, hjust=.5),
        strip.text = element_text(size=11, vjust=.7),
        plot.margin =       unit(c(1, 1, 0.5, 0.5), "lines"),
        strip.background = element_rect(colour="white", fill="white"), 
        axis.text.y = element_text(size=12), axis.text.x = element_text(size=11), 
        axis.title.x = element_blank(), axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major =  element_line(colour = "grey90", size = 0.2))
dev.off()



##################################### 
##################################### APPENDIX: LPM Tables (Appendix)
##################################### 

sink("full_sample_tables.tex")
for(i in dim(overallmeans)[1]:1){
  mod1<-lm(as.formula(paste(overallmeans$cat[i],"~ rep + ind")), data=survey)
  mod1$se<-vcovHC(mod1)
  mod2<-lm(as.formula(paste(overallmeans$cat[i],"~ rep + ind + scale(age) + income5 + income4 + income3 + income2 + black + latino + otherrace + scale(edu) + male")), data=survey)
  mod2$se<-vcovHC(mod2)
  tabb<-apsrtable(mod1, mod2, se="robust",
                  multicolumn.align = "right",
                  caption=paste("Dependent Variable: Mentioned",overallmeans$name[i]),
                  label=paste0("tab_",overallmeans$name[i]),
                  notes=list(se.note, stars.note, "Partisan leaners are coded as partisans. Age is standardized. Education has six levels and is standardized."),
                  coef.names=c("Intercept","Republican","Independent","Age","Income Quintile 5","Income Quintile 4","Income Quintile 3","Income Quintile 2","Black","Latino","Other Non-white","Education","Male")
  )
  print(tabb)
}
sink()



               